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Using the parallel tempering algorithm and GPU accelerated techniques, we have performed 
large-scale Monte Carlo simulations of the Ising model on a square lattice with antiferromagnetic 
(repulsive) nearest-neighbor(NN) and next-nearest-neighbor(NNN) interactions of the same strength 
and subject to a uniform magnetic field. Both transitions from the (2 x 1) and row-shifted (2 x 2) 
ordered phases to the paramagnetic phase are continuous. From our data analysis, reentrance 
behavior of the (2 x 1) critical line and a bicritical point which separates the two ordered phases at 
T=0 are confirmed. Based on the critical exponents we obtained along the phase boundary, Suzuki's 
weak universality seems to hold. 
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I. INTRODUCTION 

For the nearest-neighbor (NN) Ising antiferromagnet 
on the square lattice in a uniform magnetic field, the low 
temperature ordered phase is separated from the para- 
magnetic phase by a simple, 2nd order phase bound- 
ary. (Within the context of the lattice gas model 
this system could be described as having repulsive NN- 
coupling and forming a c(2 x 2) ordered state.) With the 
addition of repulsive (antiferromagnetic) next-nearest- 
neighbor (NNN) interactions the situation becomes more 
complicated. Early Monte Carlo simulations suggested 
that a single, super-antiferromagnetic, or (2x1), phase 
existed, separated from the paramagnetic phase by a sin- 
gle phase boundary d, 0, Q. A degenerate, row-shifted 
(2 X 2) state was also predicted at zero temperature. (See 
Fig.[l]for a schematic representation of these states.) On 
the other hand, symmetry arguments based on Landau 
theory' 1] predict the order-disorder transitions of (2 x 1) 
and (2 x 2) structures belong to XY model with cubic 
anisotropy. 

The original motivation of this study was to investi- 
gate the possibility of XY-like behavior of the Ising spins 
on the square lattice, since there is numerical evidence^ 
that Ising antiferromaget with attractive NNN interac- 
tion on the triangular lattice has an XY-like intermediate 
state between the low temperature ordered state and high 
temperature disordered state. In fact, the present model 
has already been studied by many authors using different 
approaches. An early Monte Carlo studyQ comprehen- 
sively showed the phase diagrams for several different in- 
teraction ratios (R) of NNN to NN interaction. But due 
to the deficiencies of computer resources at that time, 
for the R = 1 case, a disordered region was missed be- 
tween the two ordered phases which was pointed out in a 
later interfacial free energy study [fj,]. Meanwhile, trans- 
fer matrix studies 0, S] found reentrant behavior for the 
(2 X 1) transition lines. While a study using the clus- 



ter variation method[9|, |10| concluded that for a range 
of R (0.5^^1.2) the system undergoes a first order tran- 
sition, a recent Monte Carlo study[TT| using a variant 
of the Wang-Landau method[l2] focused on the R = 1 
case without external field and found the phase transi- 
tion is of second order. For external field = 4 (in the 
unit of NN interaction constant) the two ordered phases, 
namely (2x1) and row-shifted (2x2), are degenerate at 
zero temperature so it is tempting to think that the cu- 
bic anisotropy would be zero for this field and that there 
could be a Kosterlitz-Thouless transition. 

In this paper, we carefully study the location of phase 
boundaries and the critical behavior for the case R — 1. 
In Sec. In] the model and relevant methods and analysis 
techniques are reviewed. Our results are presented in 
Sec. mil along with finite size scaling analyses, and we 
summarize and conclude in Sec lIVI 

II. MODEL AND METHOD 

A. The model 

The Ising model with NNN interaction is described by 
the Hamiltonian 

71 = Jnn ^ o-iCrj + JNNN ^ a-,;crj +i7 ^ cr,;, 

(1) 

where ai,aj = ±1, Jnn and Jnnn are NN and NNN 
interaction constants, respectively, H is an external mag- 
netic field, and the sums in the first two terms run over 
indicated pairs of neighbors on a square lattice with peri- 
odic boundary conditions. Both Jnn and Jnnn are pos- 
itive (antiferromagnetic) and the ratio R ~ Jnnn/ Jnn- 
For the R = 1 case, the ground states would be the 
(2 x 1) state, also known as super-antiferromagnetic state, 
in small magnetic fields; and at higher fields it would be 
a row-shifted (2 x 2) state, which differs from the (2 x 2) 
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where N is the total number of spins and T is the sim- 
ulation temperature. In some cases, the true ordering 
susceptibility x"*", which is < (Af™^)^ >, is used to 
eliminate simulation errors resulting from < M™'^ >, 
where the order parameter is known to be zero for the 
infinite lattice. 
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FIG. 1: Schematic plots of c(2 x 2), (2 x 1) and row-shifted 
(2 X 2) ordered structures within the context of the lattice gas 
model. (In magnetic language, filled circles correspond to up 
spins and empty circles correspond to down spins.) 

state in the sense that the antiferromagnetic chains in 
the former state can slide freely without energy cost. See 
Fig. m Locally, the structure may appear to be (2 x 2), 
but for large enough lattices the equilibrium structure al- 
ways shows row shifting. As a result, such a row-shifted 
(2 X 2) state is highly degenerate, and the antiferromag- 
netic sublattice exhibits only one dimensional long range 
order. In terms of the sublattice magnetizations 

AfA-^^a„ A = 1,2,3,4 (2) 

we can define two components of the order parameter for 
the (2x1) state 

A'P = [Ml + M2 - (Ms + M4)]/4, (3) 

M'' = [Ml -I- Mi ~ (M2 + M3)]/4, (4) 

with a computationally convenient root-mean-square or- 
der parameter 

M™'^ = y^(M'^)2 -I- (M^)2. (5) 

Since M''™* would have a limiting value of ^ for the row- 
shifted (2 x 2) state and be zero for the disordered state, 
it can also be used as an order parameter for the row- 
shifted (2 X 2) state. 

Other observables, such as the finite lattice ordering 
susceptibility x a-nd fourth-order cumulant (7, are defined 
in terms of the order parameter M™^ as 

N 

— (Af™'')2 > - < Af™'' >2] (6) 

3<(M™^)2>2 



B. Simulation methods 

For small lattice sizes, Wang-Landau sampling [l3| was 
used to obtain a quick overview of the thermodynamic 
behavior of our model. A two-dimensional random walk 
in energy and magnetization space was performed so that 
the density of states g{E, M) could be used to determine 
all thermodynamic quantities (derived from the parti- 
tion function) for any value of temperatures and exter- 
nal field. Consequently, "freezing" problems are avoided 
at extremely low temperatures. This allowed us to deter- 
mine the "interesting" regions of field-temperature space; 
however, it quickly became apparent that, because of 
subtle finite size effects, quite large lattices would be 
needed. Unfortunately, as L increases, the number of 
entries of histogram explodes as L** and it proved to be 
more efficient to use parallel tempering instead. 

Since a large portion of interesting phase boundary 
is at relatively low temperatures and many local energy 
minima exist which makes the relaxation time rather 
long, the parallel tempering method [3, [l3| is a good 
choice for simulating our model. The basic idea is to 
expand the low temperature phase space by introducing 
configurations from the high temperatures. So, many 
replicas at different temperatures are simulated simulta- 
neously, and after every fixed number of Monte Carlo 
steps, a swap trial is performed with a Metropolis-like 
probability which satisfies the detailed balance condition. 
The transition probability from a configuration simu- 
lated at temperature f3m to a configuration Xn simulated 
at temperature /?„ would be 

W{X^,P„,\Xn,P„) =rmn[l,exp{-A)], (8) 

A = {f3n- (3^){n,n-Hn). (9) 

We chose the temperatures for the replicas to be in a 
geometric progression [l^, which would make acceptance 
rates relatively constant among neighboring temperature 
pairs, and the total number of temperatures was chosen 
to make the average acceptance rate above 20%. 

The multiplicative, congruential random number gen- 
erator RANECU was used [T^ . [l7j , and some results were 
also obtained using the Mersenne Twister fis'l for com- 
parison. No difference was observed to within the error 
bars. 

Typically, data from 10^ to 10^ MCS were kept for each 
run and 3 to 6 independent runs are taken to calculate 
standard statistical error bars. For parallel termpering, 
the swap trial was attempted after every MCS. In all the 
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plots of data and analysis shown in following sections, if 
error bars are not shown they are always smaller than 
the size of the symbols. 

In general, such replica exchange not only applies to 
the temperature set, but also can apply to any other sets 
of fields, such as the external magnetic field. Following 
the same argument, one can get the transition probability 
from {Xjn,H,n} to {X„,i/„} 

W{X„,,H,-^\X,-„ H„) = mm[l, expi-A)], (10) 

A = /3(if„-ff,„)(Af™-Af„). (11) 

where Mm, Mn are the uniform magnetizations of replica 
m and n, respectively. 

C. Finite-size scaling analysis 

To extract critical exponents from the data, we per- 
formed finite-size scaling analyses along the transition 
lines. Since the maximum slope of the fourth-order cu- 
mulant U follows [Toj 

(^)ma. =a'i*(l + fo'i-"), (12) 

where K = ^" , the correlation length exponent v can 
be estimated ciirectly. 

With the exponent and critical temperature Tc at 
hand, the critical exponent f3 and 7 can be extracted 
from the data collapsing of the finite-size scaling forms, 

M = L--X{tL-) (13) 



xT^LiY{tLi) (14) 

where t — |1 — and X and Y are universal func- 
tions whose analytical forms are not known. One can 
also estimate the exponent a from the relation of the 
peak position with lattice size for the specific heat 

Cmax = cL~ + Co (15) 

where Co is the "background" contribution. In some 
cases when the appropriate paths, i.e. which are perpen- 
dicular to the phase boundary, are ones of constant tem- 
peratures, then the critical behavior would be expressed 
in terms of reduced field h = \1 — and all the above 
analysis still applies. 

D. GPU acceleration 

General purpose computing on graphics processing 
units (GPU) attracts steadily increasing interest in sim- 
ulational physics [20I [2ll . [23 |. since the computational 
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FIG. 2: The phase diagram for the Ising square lattice with 
antiferromagnetic nearest- and next-nearest-neighbor interac- 
tions in a magnetic field for i? = 1. Open circles and pluses 
denote simulation results. The solid lines are second order 
transition lines, while the dashed line indicates the short range 
ordering line. 



power of recent GPU exceeds that of a central process- 
ing unit(CPU) by orders of magnitude. The advantage 
continues to grow as the performance of GPU's doubles 
every year. Recently, a GPU accelerated Monte Carlo 
simulation of Ising models [2^ was performed. Compared 
to traditional CPU calculations, the speedup was about 
60 times. 

The idea behind the implementation in Ref [l^l can be 

easily extended to our model, and the parallel tempering 
algorithm is naturally realized. Initially, all the replica 
are loaded to the global memory of the GPU. For each 
replica, the entire lattice is divided into four sublattices, 
then spins in the same sublattice can be updated simul- 
taneously by the GPU using a Metropolis scheme, and 
the swap of configurations of replica pairs can also be 
achieved in parallel. 

On the GeForce GTX285 graphics unit, our code runs 
about 10 times faster than it does on the 32 CPUs of 
IBM p655 cluster using Message Passing Interface(MPI) 
for parallel computation. 



III. RESULTS AND DISCUSSION 

A. Phase diagram and short range ordering 

From the ground state analysis, with zero or low field 
the ordered state would be the superantiferromagnetic, 
or (2x1), structure. As the external field increases to 
4 < H/Jnn < 8, more spins align in the opposite field 
direction, and the ordered structure would be the row- 
shifted (2 X 2). With even stronger fields, the system 
becomes paramagnetic. In the region near H / J^n = 4 
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FIG. 3: Variation of the specific fieat C versus field H witfi 
lattice sizes L = 20, 40, 80, f 60, 200, 300 for paths of constant: 
{a.)kBT/JNN = 1.2 and {h)kBT/JNN = 1.1. 
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FIG. 4: Nearest- and next-nearest-neighbor correlation func- 
tions. The field is varied for paths of 3 different temperatures: 
ksT/JNN = 0.4,0.5 and 0.6 across the short range ordering 
line. 



a mixture of (2 x 1) and row-shifted (2 x 2) is visible. 

For finite temperatures, we found that the fourth-order 
cumulant is always a good quantity to use to locate 
the transition points, while the data for other quanti- 
ties, such as the specific heat or susceptibility, may look 
"strange" due to the effect of neighboring critical points. 
To help the reader understand the observed thermody- 
namic properties, the final phase diagram for i? = 1 is 
plotted in Fig. [2| The solid lines are the phase bound- 
aries, all of which are continuous. The dashed line inside 
the (2 X 1) ordered phase indicates a "short range or- 
dering" line, which was located from the peak position 
of the specific heat. As shown in Fig. [3l for paths of 
constant temperature ksT/ Jnn = 1-1 and 1.2, there are 
two peaks, and the one that increases with lattice size 
corresponds to the critical point. 

An indication of the complexity of the finite size be- 
havior is clearly seen in the bottom portion of Fig. [3] in 
which the small lattices actually have minima in the spe- 
cific heat for field values that eventually show phase tran- 



sitions for sufficiently large systems. The round-shaped 
size-independent peak is due to the short range order- 
ing of the (2 X 1) "clusters" of different orientation from 
the ordered background. No corresponding behavior was 
observed from the susceptibility or the fourth-order cu- 
mulant. 

In order to confirm the above argument, we also calcu- 
lated the NN and NNN pair correlation function, that is 
< aiCj >, for paths of different fields crossing this line. 
The correlation function data are plotted in Fig. ID The 
NN pair correlation decreases from zero to a minimum 
and then increases to positive values, while the NNN pair 
correlation increases monotonically from —1. 



B. Critical behavior 

The data for the specific heat and susceptibility for 
three different values of H are plotted in Fig. [5l 

Without the field, they both show very sharp peaks, 
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FIG. 5; Specific heat and susceptibility for 3 different fields 
across the phase boundary. Data are for: L=100, A; L=120, 
• ; L=160, x; L=200 A; L=300 +; L=400 *. 

and from the magnitudes of the peak values of the spe- 
cific heat, as shown in Fig[6fb), we can obtain a rather 
accurate estimate of the exponent ratio a/v — 0.357(8), 
which is obviously not zero. In Fig|Sfa), we also show 
the curve-fit for the maximum slope of ^ for H — 0, 
and extract the exponent v = 0.847(4) directly. 

Both values of a and v are quite consistent with the 
early estimates in ref and the value of a/v is differ- 
ent from ref fll'|, in which the 1/L correction term was 
assumed up to lattice size L = 160. 

The same procedure was repeated for H/ J^n — 2.5 
and 3.3, however, as shown in Fig. [H the peaks of the 
specific heat become increasingly rounded as the field 
increases, which makes it rather difficult to get a di- 
rect estimate of the exponent a. Because of this it was 
necessary to obtain data for much larger lattice sizes, a 
task that was only possible with the use of GPU com- 
puting. In fact, as the value of v increases with the 
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FIG. 6: Curve fits using the leading terms of equations [12] 
and [15] for the maximum slopes of ^ (a) and peak values of 
specific heat (b), respectively. 

field, for H/Jj\jn = 3.3, according the hyper-scaling law 
a = 2 — dv, where d = 2 is the dimension of the sys- 
tem, a should be negative. Although the curve-fit is not 
stable, given the value of a we can get a fit within error 
bars. Such continuous increasing of the exponent v up to 
values much greater than one is actually consistent with 
the findings of an early transfer-matrix study Q. 

To estimate the critical exponents (3 and 7, we per- 
formed data collapsing with a large range of lattice sizes 
for the order parameter and its susceptibility. As shown 
in Fig.[71 the data in both finite size scaling plots collapse 
very nicely onto single curves, and the ratio [3 /v and ^ /u 
agree with values of the 2D Ising universality class within 
error bars. 

Hence, although the individual exponents are non- 
universal, Suzuki's weak universality holds quite well. 
Another data collapsing along the path of constant 
H/ Jnn ~ 6 across the phase boundary of the row-shifted 
2x2 state is shown in Fig. [8] The quality of the data 
collapsing is also excellent, and again, the exponents are 
non-universal. The estimate for (j/v is, a, bit low but ^/v 
agrees well with prediction of weak universality. 

In Table. HI the critical points and exponents a, (3, v 
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FIG. 7: Finite size scaling data collapsing along paths of con- 
stant H/ Jnn = and 2.5 for root-mean-square order param- 
eter and its ordering susceptibility, respectively, = 1 1 — 
and t = {l- ^). Data are for: L=80, o; L=100, A; L=120, 
• ; L=160, x; L=200 A; L=300, +; L=400, *. 
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FIG. 8: Data collapsing along the path of constant H/Jnn = 
6 for root-mean-square order parameter and its ordering sus- 
ceptibility, respectively. Data are for: L=80, o; L=100, A; 
L=120, •; L=160, x; L=200 A; L=300, +; L=400, *. 



order path Tc or He a 
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H=0 2.0820(4) 0.302(7) 0.103(3) 1.482(7) 0.847(4) 
2 X 1 H=2.5 1.6852(3) 0.104(19) 0.118(3) 1.657(6) 0.947(7) 
H=3.3 1.3335(6) 0.130(5) 1.930(6) 1.102(8) 



2x2* H=6 0.7293(7) 
T=0.5 6.848(5) 



0.110(5) 2.072(6) 1.176(9) 
0.126(4) 1.775(5) 1.02(2) 



TABLE I: Values of critical point temperatures or magnetic 
fields and corresponding critical exponents for several paths 
of constant temperature or field across the phase boundary of 
the (2 X 1) and *row-shifted (2 x 2) ordered phases. 



and 7 for several typical paths of constant H or T across 
the phase boundary are listed. 



C. Reentrance behavior 

Close to the region between the two ordered phases the 
correlation length exponent u turns out to be quite large, 

and correspondingly the location of the critical points 
becomes very difficult to determine. In addition, the 
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FIG. 9: Fourth order cumulant U versus field H along 
the path of constant ksT/JNN = 0.7 for lattice size L — 
32,64,128,256 . 



specific heat curves look " strange" , see Fig. [U because 
the exponent a/v would have a negative value with large 
magnitude, much larger lattice size is needed to approach 
to the limiting peak value. Since Suzuki's weak univer- 
sality seems to hold along the transition line, we fixed 
the values of P/i^ = 0.125 and ^jv ^ 1.75 for the data 
collapsing analysis to get a better estimate of the critical 
point and the exponent v. 

As shown in Fig. [51 the crossing point of the fourth 
order cumulant curves for a path of constant ksT I Jnn = 
0.7 is slightly above H/J^m = 4, and from the data 
collapsing analysis, see Fig. [TOl we obtained an estimate 
of the critical point to be He/ Jnn = 4.052(7). Hence, we 
confirm the reentrant behavior of the (2x1) transition 
line. 

For the paths of constant ksT/ J^n — 0.45, the curves 
of the fourth order cumulant shows two crossing points 
and the finite size effect is quite obvious. See Fig. [TT] 

For the larger lattice size, the two crossing points move 
towards lower fields but they do not approach each other. 
Thus, a region of disorder remains between the two differ- 
ent ordered states, even down to quite low temperatures. 
(If however, small lattices are used with insufficient data 
precision, it looks as though the curves for different lat- 
tice sizes coincide. Such behavior would indicate, erro- 
neously, the presence of an XY-like region.) In Fig. [T^l 
we show data collapsing fits along the path of constant 
kBT/ Jnn = 0.5 (and excellent data collapsing is also 
found along the path of constant ksT/JNN — 0.3), which 
confirms that there is a disordered region in between the 
two ordered states. 

We thus conclude that there is no XY-like region and 
that the two phase boundaries probably only come to- 
gether at a bicritical point at T = 0, although we can- 
not exclude the possibility of a bicritical point at very 
low, but non-zero, temperature. However, the data do 
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FIG. 10: Finite size scaling data collapsing along the path 
of constant ksT/ Jnn = 0.7 for root-mean-square order pa- 
rameter and its susceptibility, respectively, ft' = 1 1 — | and 
ft = (1 - ■^). Data are for: L=32, <; L=64, L=128, >; 
L=256, <. 
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FIG. 11: Fourth order cumulant U versus field H along 
the path of constant ksT/JNN = 0.45 for lattice size L = 
64,128,256. 
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The variation of the critical exponents is consistent with 
the changing magnetic field producing different effective 
anisotropies which, in turn, is expected to yield non- 
universal behavior 23] . Due to the large values of v near 
the bicritical point (and correspondingly strongly nega- 
tive values of a), we consider it also extremely unlikely 
that tricritical points could be found along these transi- 
tion lines as T becomes small. 
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FIG. 12: Finite size scaling data collapsing along paths of 
constant ksT/JNN = 0.5 for root-mean-square order param- 
eter and its susceptibility, respectively. Data are for: L=80, 



not yield any hint of such a bicritical point; but the 
lack of data points for T < 0.2 in Fig. 2 precludes us 
from making a definitive statement about this issue. 
(Moreover, the reentrant behavior of the (2x1) phase 
makes it very difficult to study the approach to T = 0.) 



IV. CONCLUSION 

We have carried out large-scale Monte Carlo simula- 
tions for the square-lattice Ising model with repulsive 
(antiferromagnetic) nearest- and next-nearest-neighbor 
interactions. From the finite size scaling analysis, both 
transitions from (2x1) and row-shifted (2 x 2) ordered 
states to disordered states turn out to be continuous and 
non-universal. The reentrance behavior of the (2 x 1) 
transition line is confirmed, and the proximity to the 
transition line to the (2 x 2) state makes it difficult 
to untangle the low temperature behavior unless quite 
large lattices are used. It was only possible to obtain the 
precise, large lattice data needed though the use of GPU 
computing. No evidence for XY-like behavior is found, 
and we conclude that there is probably a zero temper- 
ature bicritical point. Although the exponent v varies 
along the transition line, the exponent ratio and j /v 
seem to agree with that of the 2D Ising universality class. 
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